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ABSTRACT 

Chaotic flow is studied in a series of numerical magnetohydrodynamical simula- 
tions that use the shearing box formalism. This mimics important features of local 
accretion disk dynamics. The magnetorotational instability gives rise to flow turbu- 
lence, and quantitative chaos parameters, such as the largest Lyapunov exponent, can 
be measured. Linear growth rates appear in these exponents even when the flow is 
fully turbulent. The extreme sensitivity to initial conditions associated with chaotic 
flows has practical implications, the most important of which is that hundreds of or- 
bital times are needed to extract a meaningful average for the stress. If the evolution 
time in a disk is less than this, the classical a formalism will break down. 
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1 INTRODUCTION 

Astrophysical accretion disks are able to evolve because an- 
gular momentum is extracted from fluid elements and trans- 
ported outward. This is effected by the presence of non- 
vanishing radial-azimuthal components of the Maxwell and 
Reynolds stress tensors, produced by magnetohydrodynamic 
(MHD) turb ulence driven by the mag netorotational insta- 
' bility (MRI) jBalbus fc Hawlev l ll998L As there is no ana- 
lytic theory of MHD turbulence at hand (nor is there one in 
sight), large-scale numerical simulation has been the main 
avenue of progress toward understanding its properties. 

Many numerical st udies make use of the local "s hear- 
ing box" approximation jHawlev. Gammie fe Balbus Il995t) . 
The shearing box is an invaluable tool for studying local flow 
dynamics in detail, and for resolving turbulent flow with the 
largest possible dynamical range. MHD turbulence, like all 
true turbulence, should be chaotic, as quantified formally 
by a measured positive Lyapunov exponent. What is less 
clear, but of great astrophysical significance, is whether long 
term flow averages are even well-defined. To put the ques- 
tion most starkly, imagine macroscopically identical disks, 
with fluid perturbations that vary by a tiny amount. The 
fine scale description of their internal turbulence will surely 
differ, but will quantities such as the stress tensor compo- 
nents converge to the same values in the long term? This 
question g oes directly to the heart of the standard a disk 
formalism JShakura fc Sunvaev lll97l which assumes that 
disk transport may be described by a spatially constant or 
slowly varying a parameter. (The quantity a is defined as 
the ratio of radial-azimuthal component of the stress tensor 



Tnq> to the gas, or gas plus radiation, pressure.) How well 
supported is this assumption? 

This paper begins to study these complex issues by ex- 
amining the chaotic nature of the MHD turbulence. The pre- 
sentation is organized as follows. In §2 we briefly review the 
shearing box. In §3 we present the results of experiments 
designed to reveal how the measured properties of turbu- 
lent flow are related to both computational and physical in- 
put parameters. In §4, we carry out a series of experiments 
that demonstrate qualitatively that the MHD turbulence is 
chaotic, and in §5 this is quantified by computing the Lya- 
punov exponents for a set of simulations. Our conclusions 
are summarized in §6. 



2 SHEARING BOX SYSTEM 

The shearing box syste m, described in 

IHawlev. Gammie fc Balbus I dl995l) . is designed to rep- 
resent a very local section of an accretion disk, viewed in 
corotating coordinates with angular frequency Q. Starting 
with the full set of dynamical equations in cylindrical 
coordinates, the equations are locally expanded about a 
fiducial cylindrical radius R, with (dx, dy, dz) corresponding 
to cylindrical coordinates (dR,Rd(j>,dz). The computational 
domain is then a Cartesian box, but with the rotational 
inertial forces (Coriolis and centrifugal) retained. The 
centrifugal force nearly balances gravity, leaving a remnant 
tidal force linear in x. All other forces are directly retained. 
For this system, the ideal MHD equations of motion 
become, 
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retain their usual forms. The internal energy per particle e 
is defined by 



P = pe( 7 -1). 



(5) 



For an isothermal gas, the energy equation is replaced with 

Pp~ X — constant. (6) 

The variables in these equations have their usual meanings. 
We shall also introduce the constant q, which is a param- 
eter describing the local radial dependence of the angular 
frequency, i.e., q — — dlntt/d\n R. For a Keplerian angular 
momentum distribution, q — 3/2. In this study we have ig- 
nored vertical gravity, and have thus dropped the Q 2 zz term 
from equation (0. 

These equations are solved using the time-explicit, 
operator-split finite differe ncing algorithm o f the ZEUS 
code f or hydrodynamics l|Stone fc Norman | | l992cj) an d 
MHD JStone fc Norman I Il992bl: lHawlev fc Stone I Il99sh . 
The shearing box computational domain extends in x, y, 
and z respectively from —L x /2 to +L x /2, to L y , to L z . 
The boundary conditions are periodic in both the y and z 
directions, and "shearing periodic" in x, meaning we equate 



f{x,y,z) = f{x + L x ,y- qflL x t,z). 



(7) 



The azimuthal velocity has an additional correction to offset 
the angular velocity difference between the inner and outer 
radial boundaries: 



v y (x,y, z) = v y (x + L x ,y - qttL x t, z) + qVLL x 



(8) 



The initial state is one of uniform density and pres- 
sure, and the velocity profile is v = —qQxy. In this paper, 
our initial magnetic field configurations are either a uniform 
toroidal magnetic field, or a vertical magnetic field vary- 
ing sinusoidally in the radial x direction. The magnetic field 
strength is described by the standard plasma f3 parameter 



8ttP 



(9) 



the ratio of the gas to magnetic pressures. The MRI is trig- 
gered by seeding the initial equilibrium state with small 
pressure and velocity fluctuations. 



3 SATURATION AMPLITUDE IN THE 
SHEARING BOX: A SURVEY 

One of the goals of shearing box simulations has been to 
try to understand what determines the sustained saturation 
amplitude of the MRI-driven turbulence. In accretion disk 
applications, the focus is often on the T r ^, component of the 
combined Reynolds and Maxwell stress tensor, 



-B x B v /4tt + pv x 5v y 



(10) 



Here Sv y is the azimuthal velocity fluctuation, obtained by 
subtracting the equilibrium shear flow. This can be directly 
measured in shearing box simulations. The volume-averaged 
stress is usually normalized to the averaged gas pressu re, 
a quantity equivalent to the flShakura fc Sunvaev Ill973l) a 
parameter. 

Even with a system as simple as the shearing box, there 
are many possible factors, both computational as well as 
physical, that could in principle govern the saturation am- 
plitude. These include the initial field strength and geome- 
try, box size and aspect ratio, value of f2, background gas 
pressure, and numerical resolution. Many of these are physi- 
cally significant only in relation to others. For example, fi is 
an arbitrary time scale in the system. The shearing box sys- 
tem is independent of O so long as q remains constant and 
the physical box size changes proportionately. The formal 
scaling invariance is that if 



(11) 



where c is an arbitrary constant, then the equations remain 
unchanged if 



(x,t) -> (x/c,t/c). 



(12) 



We will make use of this fundamental property in §4 when 
we explore chaotic behavior. 

For the cases of an initial uniform field oriented along 
either the y or z axes, lHawlev. Gammie fc Balbus I dl995l) 
found that the saturation amplitude of the turbulence de- 
pended upon both the box size and initial field strength. 
When the field had a vanishing mean value over the compu- 
tational box volume, the final outcome was much less sen- 
sitive to the initial field energy. Here, we present a series 
of zero mean field simulations designed to test how the box 
size affects the overall level of the turbulence. The baseline 
calculation has q = 1.5 corresponding to a Keplerian back- 
ground, and a box size of L x = 1, L v = 2n, and L z — 1, 
set so that L z = c 2 /Q where c 2 — P/p. The actual val- 
ues used are P = 10 -6 and Q = 0.001. The gas is adia- 
batic with 7 = 5/3. The initial field is vertical and varies 
sinusoidally in the x direction, B z oc sm{x / L x ). The field 
strength set so that the volume average magnetic energy 
corresponds to /3 = 800. The maximum vertical field has 
(3 — 400, which gives a fastest growing wavelength for the 
MRI of A max = 8tt/15 1/2 (v A /ty = 0.46. The grid resolution 
is 32x 64x 32 zones in (a;, y, z). This setup corresponds to the 
so-cal l ed "fi ducia l box" used in lHjwlev.Gammi e fc Balbus I 
Jl995lll996l) and lBrandenburg et al. I (Il995f) . 

In addition to the fiducial run, we computed four ad- 
ditional simulations with different box dimensions. In three 
simulations we doubled each of the box dimensions in turn, 
and in the final simulation of this set we quadrupled the ver- 
tical size. To maintain a constant resolution number of grid 
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Figure 1. Magnetic energy histories for an ensemble of shearing 
box runs with different box sizes. Magnetic energy is expressed in 
terms of the ratio of the magnetic pressure to initial gas pressure, 
fi~ 1 . The zero point in time for each run is offset by 30 orbits for 
ease of comparison. Each curve is labeled by the dimension that 
is increased. The fiducial box is the rightmost plot. 



Figure 2. Magnetic energy histories for an ensemble of shearing 
box runs with different vertical box sizes. Magnetic energy is ex- 
pressed in terms of the ratio of the magnetic pressure to initial 
gas pressure, /3~ . The zero point in time for each run is offset 
by 30 orbits for ease of comparison. Each curve is labeled by the 
value of L z . 



zones per unit length, the number of grid zones used was in- 
creased appropriately. Figure shows the magnetic energy 
per unit volume evolving over 30 orbits for each run in this 
series of models. While one must be cautious extrapolating 
from this limited time baseline, the results generally suggest 
that increasing the vertical dimension is the most important 
effect over the long term. The temporal variations in mag- 
netic energy within a given run are, however, quite large, 
often well in excess of the running average. 

To isolate the effect of the vertical size, we carried out 
a second set of simulations in which L z is varied, from 
L z = 1/4 to L z = 4, and all other factors are held fixed, 
including the grid zone size Az. The initial magnetic field 
is as described above. Here we used an isothermal equation 
of state, 7 = 1.0 which keeps the box from heating due to 
dissipation of the turbulence. Figure [5] shows the time evo- 
lution of the magnetic energy and even over a limited time 
baseline clearly reveals a general trend toward increasing 
magnetic energy with increasing z box size, again with sig- 
nificant variation. The most striking conclusion is that there 
is a lower bound to the vertical size required for amplifica- 
tion to occur: the field energy declines dramatically in the 
L z — 1/4 model, and even the initial linear growth rate is 
reduced. The fastest growing vertical wavelengths for initial 
field strengths with /3 < 1350 no longer fit inside the box, 
and evidently losses win over the reduced power input at the 
top of the turbulent cascade. 

For box sizes L z = 1/2 and larger the turbulence is 
sustained over the evolution period of 30 orbits. However, 
other trends are visible. First, the amplitude of the initial 
spike of linear growth saturation decreases with increasing 
L z . Second, the amplitude of the magnetic energy variations 
in the turbulent state is reduced with increasing box size. 

These history plots demonstrate that it is not a sim- 
ple matter to quantify the saturation amplitude of the tur- 
bulence in a shearing box, or to determine unambiguously 
the influence of physical parameters in setting that ampli- 



tude. It is not obvious, in general, over what length of time 
one should average to obtain a characteristic amplitude; nor 
is it even obvious that such a baseline must always exist. 
Large fluctuations, even in volume- averaged data, are the 
rule. This is the overarching and defining characteristic of 
the MHD turbulence: it is chaotic. 



4 ONSET OF CHAOS 

It is clear that changing the physical parameters of a shear- 
ing box run can result in different levels of turbulence. Here 
we show that almost exactly the same physical parameters, 
but with tiny variations in numerical value, leads to macro- 
scopically different levels of turbulence. The sensitivity of 
the system to infinitesimal perturbations is a key feature of 
chaos. 

The first set of simulations consists of three shearing 
boxes with different SI values, but which have been rescaled 
so as to be physically equivalent. All three have q — 1.5, 
/3 = 800, 7 = 1.0, B z oc sm(x/L x ), and 32 x 64 x 32 grid 
zones. One run uses the standard box size of 1 x 2n x 1 
with n = 0.001. The other two have Q = 0.0031 and 
SI — 0.000507 with the physical box sizes scaled accord- 
ingly. 1 An identical sinusoidal velocity perturbation was ap- 
plied to the initial state in each box. In principle, since 
the length to time ratio is preserved between these differ- 
ent shearing box systems, the resulting evolution should be 
identical. In the computational experiment, however, these 
scaled systems evolved very differently. 

Figure |3] shows the volume averaged magnetic energy 
density histories for the three runs. The initial linear growth 



1 Care must be taken so that the rescaling does not correspond 
to a simple shift by an integer number of bits in the machine 
registers. 
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Figure 3. Time evolution of the magnetic energy for a series 
of simulations with different Q values, with the box dimensions 
adjusted so that each run is formally equivalent. Despite this, 
the magnetic energies diverge after the initial linear growth and 
saturation. 



Figure 4. Running time averages of a for the three simulations 
with different Q values. The averages are computed beginning 
with orbit 45, well after turbulence is established. This shows 
that the average a values only slowly converge to comparable 
time-averages even after 120 orbits. 



and saturation phase is indistinguishable. Significant devi- 
ations appear after orbit 8, and at orbit 20 the energies 
are widely different. Since these simulations correspond to 
identical physical conditions, one might expect that, at the 
very least, time-averaged values would coincide. Figure 2] 
illustrates how the running time averages of the volume- 
averaged stress parameter a for each run evolves with time. 
The running time average is defined 

<«(«)> = ^ J\dt, (13) 

and the expectation is that this should converge to a steady 
value characteristic of the turbulence. The question is, over 
what time interval should this happen? As figure 0] shows, 
the fluctuations in the stress are so large that the running 
time averages (a(t)) do not coincide. Since a very large ex- 
cursion in magnetic energy is visible in figure|^]at 20 orbits, 
we computed the running average from 45 to 120 orbits. But 
over this interval a unique characteristic a does not emerge 
even for a single simulation. 

The next experiment further demonstrates the chaotic 
nature of the turbulence by taking the data halfway through 
an extended run and randomly perturbing the flow veloci- 
ties with a gaussian rms value of 0.001%. The subsequent 
evolution is then compared with the extended run. The pa- 
rameters for this test were Q = 0.0005, q = 1.5, /3 = 400, 
7 = 1.0, B z oc sin(x/L x ), grid size of 32 x 64 x 32, and 
an appropriately scaled box equivalent to the standard box. 
The results are shown in figure |S] The arrow indicates the 
time when the perturbations were added, and the perturbed 
simulation was then run on to 60 orbits in time. The plot of 
magnetic energy reveals a visible divergence by orbit 30, and 
a substantial difference by orbit 35. This is the denning fea- 
ture of chaos: substantially different macroscopic behavior 
created by infinitesimal perturbations. 

To compare the behavior of models with different ini- 
tial field topologies, we perform an experiment with an ini- 
tial uniform toroidal magnetic field for two boxes, one with 
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Figure 5. The evolution of a baseline simulation (solid line) and 
a comparison simulation to which small velocity perturbations 
are added (dashed line). Data from the baseline simulation are 
perturbed at the 0.001% level at the point in time indicated by 
the arrow. The resulting magnetic energy histories diverge visibly 
by 30 orbits, and are markedly different by 35 orbits. 

fi = 0.001 the other with fi = 0.0031. Other physical pa- 
rameters were the same as in the vertical field experiments. 
FigureHJdisplays the volume- averaged magnetic energy den- 
sity histories for these runs. Again, the initial linear growth 
phases of the two runs were identical; diverging behavior is 
a nonlinear phenomenon. The amplitude of the variations 
is not as great as in the vertical field case. Nevertheless, as 
figure |7| shows, the running time averages of a required 100 
orbits to converge. 

As a contrast to the chaotic behavior of MHD turbu- 
lence, we also ran several purely hydrodynamic shearing box 
systems. Hydrodynamic shearing boxes are linearly stable 
when q < 2 (the Rayleigh criterion), and nonlinearly sta- 
ble when q < 2 — e, where e is a number on order 0.01 
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Figure 6. Volume averaged magnetic energy as a function of time 
for two different f2 simulations with an initial uniform toroidal (y) 
magnetic field. Magnetic energy is expressed in terms of the ratio 
of the magnetic pressure to initial gas pressure, 
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Figure 7. Running time averages of a for two toroidal field sim- 
ulations with different Q values. The averages are computed be- 
ginning with orbit 45, well after turbulence is established. The 
time averages approach different values by 120 orbits. 



jHawlev. Balbus. fc Winters Ill999t) . All of the q = 1.5 hy- 
drodynamical simulations gave similar results. The volume- 
averaged kinetic energies die out in an identical fashion for 
equivalent problems that use different Q values. Small fluc- 
tuations exist in the form of waves, and are typical of de- 
caying hydrodynamic turbulence. These waves were repro- 
duced almost identically from simulation to simulation — the 
measured differences were consistent with round-off error. 
Additional experiments with different q parameters (up to 
q — 1.95) were carried out, and all of these stable hydro- 
dynamical experiments remained identical in their macro- 
scopic properties. The numerical differences were small, not 
increasing, and consistent with machine precision. These re- 
sults are consistent with the conclusion that the observed 
chaos is indeed a property of the MHD turbulence, not of 
the shearing box itself or the numerics. 



4.1 Lyapunov Exponents 

The experiments provide compelling, but qualitative evi- 
dence that MRI-driven MHD turbulence is chaotic. A more 
quantitative parameterization of chaotic behavior is afforded 
by the use of Lyapunov exponents. A system with N degrees 
of freedom has N Lyapunov exponents. The largest positive 
exponent represents the average rate of divergence of two 
initially close evolution paths in phase space. The evolution 
paths are traced out by monitoring the state vector of the 
system, which we define as 



Vx,Vy,V z , ■ 



(14) 



The density itself is generally of secondary importance, and 
therefore is tracked only indirectly in the Alfven velocity. 
Close evolution paths imply that there are small differences 
between the state vectors for each of the two trajectories. If 
an initially small difference diverges exponentially between 
two state vectors, the system is formally chaotic. This cor- 
resp onds to a positive Lyapunov ex ponent . 

iKurths fc Brandenburg I Jl99lT) demonstrated chaos in 
MHD solar convection simulations in three-dimensions and 
computed estimates of the Lyapunov exponents. We follow 
their procedure to perform a similar analysis for MRI-driven 
turbulence. We begin with an evolution in which MHD tur- 
bulence has fully developed. Then, a second evolution path 
is created by randomly perturbing the velocity components 
of the first system's state vector. The two state vector tra- 
jectories are then evolved for about 4 orbits, and, at regular 
time intervals the fractional state vector difference between 
the two paths was calculated: 



|Vp-V| 
|V| 



(15) 



where V p and V are the perturbed and original state vec- 
tors. Data from the last two orbits in the state vector dif- 
ference history are then fit to an exponential from which we 
estimate the value of the largest Lyapunov exponent. 

In the first experiment we employ a standard shearing 
box, namely 1 x 2ir x 1, grid size of 32 x 64 x 32, with an 
isothermal equation of state, 7 = 1.0, sinusoidal vertical 
field, B z oc sin(x / L x ) , of average magnetic energy f3 — 800, 
and a background Keplerian shear, q = 1.5. This was evolved 
until MHD turbulence had fully developed. Perturbations 
were then applied to a parallel computation and the state 
vector difference was calculated as a function of time. The 
resulting Lyapunov exponent is 0.458cj max , where w m ax = 
0.75Q, is the maximum unstable growth rate for the linear 
phase of magnetorotational instability l|Balbus fc Hawlevl 
1998). All of our experiments yield a Lyapunov exponent 
comparable to w ma x- 

That the Lyapunov exponent would be on order of the 
maximum MRI growth rate is not surprising. It is precisely 
the linearly unstable, exponentially growing MRI that is 
feeding the turbulence, and driving exponential divergence 
of the state vector. To examine this more carefully, a se- 
ries of experiments was performed on shearing boxes with 
a variety of background shear q parameters and for both 
sinusoidal vertical and toroidal initial field configurations. 
Because the maximum li near growth rate of the MRI is q/2 
(Bal bus fc Hawlev I1 998L the ensemble of simulations spans 
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Figure 8. The percent state vector difference histories for a series 
of vertical field simulations (solid lines) and toroidal field simu- 
lations (dashed lines). Each curve is labeled by the background 
shear parameter q. The exponential divergence that sets the first 
Lyapunov exponent is clearly seen. The rough trend of larger 
Lyapunov exponent with increased q value is also visible. 

an interesting range of MRI growth rates. We would ex- 
pect the Lyapunov exponent to be proportional to q as well. 
Figure [S] displays the volume averaged, percent state vector 
difference histories for these runs. The curves are labeled 
by their q values, with solid lines for the vertical field runs 
and dashed lines for the toroidal field runs. The correspond- 
ing first Lyapunov exponent values normalized by oj max are 
0.521, 0.458, 0.422 and 0.583 for the vertical field cases, in 
order of descending q value, and 0.484 and 0.644 for the 
toroidal field q — 1.5 and q — 1.3 runs. Clearly, the first 
Lyapunov exponents are positive, and of the same magni- 
tude, normalized by oj max . It should also be noted that first 
Lyapunov exponents were calculated at many points in time 
in the simulations, always with similar results. 

In summary, the Lyapunov exponents in MRI-driven 
MHD turbulence simulations are all positive, and, when nor- 
malized by oj max , they all lie within in a range of 0.4 to 
0.6. There is a slight trend of larger Lyapunov exponent for 
larger q value (non- normalized) . Previous experiments also 
foun d stronger overall levels of turbule nce with larger q val- 
ues ijHawlev. Balbus. fc Winters lll9 99h ; in a sense, a larger 
Lyapunov exponent is "more turbulent." Finally, the chaos 
parameters of the turbulence are independent of the initial 
magnetic field configuration. Initial vertical magnetic fields 
have similar Lyapunov exponent values as initial toroidal 
magnetic fields. 



5 CONCLUSIONS 

Shearing box simulations of the MRI constitute an excellent 
numerical laboratory in which to study chaos in turbulent 
flow. Their compactness and simple boundary conditions 
make them a very convenient system to study, but they also 
require care to interpret. For example, the variance of the 
the flow fluctuations, which may be of direct astrophysical 
interest because of its connection with radiative emission, is 
a function of the box size adopted (cf. figure [5J. 



In this paper, we have demonstrated the extreme sensi- 
tivity of MHD turbulence to infinitesimal deviations in the 
flow. This was done by several different methods: showing 
that invariant scaling laws fail when implemented numeri- 
cally, for both vertical and toroidal initial fields, and exter- 
nally imposing tiny perturbations on an established turbu- 
lent flow and following the growing deviations in the sub- 
sequent evolution of the original and the perturbed system. 
Estimates of the largest Lyapunov exponent in a variety of 
turbulent flows with different field geometries yielded values 
near the characteristic growth rate of the linear MRI. This is 
an indication that the linear physics of the instability plays 
an active role in defining the highly nonlinear turbulent dy- 
namics of these flows. One way this could come about would 
be if the energy input into a Kolmogorov-like cascade was 
essentially the linear MRI. 

The most important practical consequence of this 
behavior is that the nongaussian statistical properties 
of chaotic flows severely limit the extent to which a 
modeling can be used uncritically. Though Maxwell and 
viscous stress have some f ormal properties in common 
(Balb us fc Papaloizou II1999I) . the averaging procedure nec- 
essary for a semi-local treatment of the turbulence is a very 
delicate matter. A time base of hundreds of orbits is clearly 
necessary to establish a meaningful estimate of the charac- 
teristic stress. In astrophysical systems, especially those in 
transience, it may not be possible to ascribe an instanta- 
neous a value to the stress, and there may be no recourse 
other than detailed numerical modeling. 
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